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List of Symbols 


A: cross-sectional area of one side of the actuator 

a: assigned constant 

B: magnetic flux density in the gap 

Bm: magnetic flux density inside the actuator 

b: assigned constant 

c 2 : assigned constant 

F (x, i) : attractive force between the actuator and the rim (suspended 
mass) 

G p (s) : plant dynamic equation or open loop transfer function 

g,, 9i f 9 : 9 a P distance 

g 0 : equilibrium distance 

X : input current of the actuator 

I b , I 0 : equivalent bias current 

i: controlled current 

L: inductance of the actuator 

L c : inductance of the coil 

Lji inductance of the gap 

1: the length of actuator 

m: mass of the rim (suspended mass) 

N: n umb er of turns of he coil 
R: actuator resistance 
s: Laplace transform variable 
V: input voltage of the actuator 
V b : equivalent bias voltage 
v: controlled voltage 

x: controlled distance or small perturbation distance from 


equilibrium 

|i : permeability of the actuator material 

: permeability of free space 

4> : magnetic flux 
L.H.P : left half plane 

ASPS : Annular Suspension and Pointing System 



I. Introduction 


A. Problem Definition 

Frictionless electromagnetic suspension and levitation has 
attracted much attention since 1970. Applications include high- 
speed machine tool spindles, ultra-centrifuges, high vacuum pumps, 
and fly-wheels for energy storage. Methods of producing 
electromagnetic suspension and levitation include controlled DC 
electromagnets, diamagnetic materials, superconductors, hybrid 
systems, and tuned LCR circuits. A comprehensive review lecture of 
electromagnetic suspension and levitation techniques can be found 
in reference 1. 

The technique of suspension and levitation with controlled DC 
electromagnets is the most advanced and successful at this time. 
Many investigations are underway worldwide. Advanced ground 
transportation schemes, contactless bearings for ultra-high speed, 
and gyroscopes have been successfully demonstrated by many groups 
of researchers. 1 

The Annular Suspension and Pointing System (ASPS) developed by 
the Flight System division of Sperry Corporation 1 is a six-degree 
of freedom payload pointing system designed for use with the space 
shuttle. This magnetic suspension and pointing system provides 
precise controlled pointing in six-degrees of freedom, isolation of 
payload-carrier disturbances, and end mount controlled pointing. 
Those are great advantages over the traditional mechanical joints 
for space applications. More detail discussions of the magnetic 
suspension joints and mechanical joints can be found in reference 
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Figure 1 and 2 show the ASPS designed by Sperry Corporation. 
It consists of six actuators, three for vertical movements, two for 
radial movements, and one for tangential movements. By the coupling 
and decoupling matrices (figure 3) 2 , we can carefully decompose the 
command signal of each degree of freedom to each actuator 
individually. In other words, the coupling and decoupling matrices 
change the six— degree of freedom ASPS control system to six single- 
degree of freedom ALPS control systems. Hence, we can design each 
control loop separately. 


figure 1. 





figure 5^. 


VERNIER ACTUATOR FORCES 
AND MOMENTS 



F x “ Fy cos 450 + F V cos 450 
Fy “ Fy SIN 45° - F v SIN 45° - F w 
F Z * F A +F B +F C 

T0 =■ F b Rq SIN 60° - F c Rq SIN 60° + Fy Zq SIN 45° - F v Zq SIN 45° - F w Z 0 
JO « -F a Rq ♦ F b Rq COS 60° + F C Rq COS 60° - Fy Z 0 COS 45° - F v Z 0 COS 45° 
Ti p = Fyy Rq 


d 


VERNIER ACTUATOR DECOUPLING 




FIGURE J> 


1 





















B. Project Objectives 


(i) Understanding the existing ASPS in the lab. 

(ii) Model the dynamics of a single DC controlled ALPS 

actuator as accurately as possible. 

(iii) Re-design a controller for the single degree of freedom 
ALPS control system to achieve the highest stiffness as possible. 
{Highest stiffness will have the lowest motion in response to 
external forces} 

c. Summary 

In this design, we first analyzed the assumed model of the 
single degree ASPS bearing actuator, and obtained the plant 
dynamics eguations. By linearizing the plant dynamics eguations, we 
designed the cascade and feedback compensators such that a stable 
and satisfied result was obtained. The specified feedback 
compensator was computer simulated with the nonlinearized plant 
dynamics equations. The results indicated that an unstable output 
occurred. In other words, the designed feedback compensator is 
fail. The failure of the design is due to the Taylor's series 
expansion does not converge. 
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n. Modelling of the single degree ASPS bearing actuator 

a. Assumptions of the model 

The single degree ASPS bearing actuator consists of two pairs 
of magnetic coil elements, mounted in opposition, to control the 
rim (suspended mass) along a single axis. Figure 4 shows the 
configuration of the actuator. The magnetic coil elements have 
current biasing superimposed by a controlled DC voltage source to 
produce a force to suspend and point the rim. For fine pointing 
application, the gap distances between the rim and actuators are 

kept to a predetermined value(<0. 3") . 

figure 4. 


•ti 



ii 



We derived the equation for the single ALPS bearing actuator 

with the following assumptions: 

(i) The force of attraction betweeri magnetized bodies is given 

by F = (l/2u) * B 2 * area. 

(ii) The magnetic flux density is uniform between the gap, or 

gap distance << size of actuator. 

(iii) The rim is a perfect conductor. That is, it does not 
support any magnetic field strength (H) inside the rim. 

(iv) There is no coupling effect between each coil elements. 

(v) The controlled electromagnet behaves linearly, and there 

is no loss or fringing of magnetic field. 

(vi) The rim has a plane area over the magnetic coil 

assemblies . 

B. Derivation of ASPS D ynamics 
a) B-field in the air-gap. 

From magnetic circuit theory, the total current linked by the 
path of an N-turn coil is given by 


? 


Hili 


= N i 


( 4 ) 


Assume 

(i) the magnetic material is approximated by B m = + 

B c where B 0 is a constant. 

(ii) gap distant at path (1) = gap distance at path (2) 
ie. g, = g 2 =g 

(iii) from boundary condition, B inside the coil * B in 
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the gap. 


Therefore, equation (1) becomes, 


# 1^1 + H 29 2 + Hm 1 = Nl 


- 2Hg + 



= Ni 


(because the 
rim is a perfect 
conductor) 


— 2 Hg + —1H = Ni + 

V-m 


V-m 


As the term B 0 (1/mJ can be equivalently assigned to 
a magnetomotive force (mmf) (B 0 //x m )l 


So 


2 Hg + 


= Ni + NI a 
Vm 


Thus the magnetic flux density at the gap is given by 


B = \i a H = 


* Jq) 
2 g + -^1 


As for the particular material we used for actuator, H m »V 0 ' 
Hence 


B « 


+ i 0 ) 

2g 


( 2 ) 
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b) Relationship between F, I and g 
By considering the stored magnetic energy, Bohr 5 and Hayt 6 were 
able to relate the magnetic attraction force to the magnetic flux 
density and cross-sectional, 

F = — =— B 2 * area 

* 2 \i 0 


The geometry we used is similar to Humphris 7 and Groom 8 , figure 
4, which have two electromagnets positioned opposite the rim. This 
kind of configuration is more linear if we separate the magnetic 
flux density into the controlled and bias components. 7 

Assume no coupling effect between the two actuators, by 

equation (2) : 

B i = 2g[ 2 g 2 

Let H t = N 2 = N, I, = I 2 = I., i. = -i 2 - i = controlled current. 
Therefore, the total force acting on the rim is given by 
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F = F 1 - F 2 = <Bi - B\) = + B 2 ) (B 2 - B 2 ) 


. - M<r^ r £o — ] 

Consider B x B 2 2 - x gr 0 + x 

uJV (Sr n + x) (J p + i) - ig 0 - x) (J c - i) 1 

2 " x 2 


^ [gr 0 (I 0 + i - l 0 + i) + *(*o + i+ -To - i) 1 


2 (g 0 2 - * 2 > 


u_f/ . |i«w(sr 0 i + I#*) 

p£L (2 g a i * 21^) = Z+—2. 

2 ( gl - X 2 ) 9o - * 2 


Similiarly 


B + B = * J ° ~ i) » 

' S 2 2 L or - x cr. + X 


Mo* 


< 7 o 2 - * 2 


(g 0 J 0 + xi) 


2 g D - x g 0 

= _A ^<9^ * J o*> j [ M^£o£o_Li^_ 


Thus F = — (B x + S 2 ) (Si - B 2 ) » 2 2 

Mo Po 9o ~ x 


gl - x 2 


^.^ 2 A{g 0 i + JpX) (gr 0 J 0 + xi) 

(g c 2 -x 2 ) 2 


By the Taylor's Series Expansion at the equilibrium point (x e/ i 0 ) , 
we get 
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F(x,i) - F(x 0 , i 0 ) + (x - x a , i - i 0 )*VF| (Xo> Ig) 


a S dF = VcX 2A _ [g (g I + xi) + X(gj + J *x) 1 

32 (grj -X 2 ) 2 

M = — [( gl - x 2 ) 2 [I 0 (g 0 I 0 + xi) + i (g 0 i + Ijx)l 

dx (gl - X 2 ) 4 

- 2 (gl - x 2 ) ( ~2x) (g 0 I 0 + xi) (gr 0 i + Ijc) ] 

at equilibrium point ( x Q , i D ) = (0, 0), 

F (x 0 , i a ) = 0 

dF\ = 

dx ^ x °' *** q 3 


3F. _ 

ai IU e . i 0 > “ 2 


So 


F ( x, i) - X * |£| (jfo . i o) + i * -^1 


3F| 


5i ' x °' 


{x ot i c ) 


= < 




srl 


gl 


(3) 


c) V-I relationship of the actuator. 
Recall equation (2) : 

yJJU + J p ) 


B = 


2sr 
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Therefore, the 


By definition, *=NBA and L-(d$/di) 5 . 
inductance of the gap is 





2g 


and the inductance of the whole circuit is 


L = L, + I* + L c = 


N 2 A\l Q + N 2 A\l Q + L 
2^ 2g 2 


N 2 A\i 0 g 0 

gl - x 2 


+ L r 


(4) 


(by 4) 


By Kirchhoff's voltage law, we have 


V = Ri + (Li) 
dt 


„ . - di . . dl ^ dx 

= Ri + L-^r + j-*-j r 

dt dx dt 


We previously separated the voltage, current and gap distance 
into the bias components and controlled components. 

That is Let V = V b + v 

i - I k + i 

and x = g 0 + x 

Therefore, 
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. . . _ di . dl . dx T dl dx 

* V =■ R(I b * 1) + L-^ * 1 -£ c , Tt b ~dZ ~di 


As, at equilibrium position, V b = RI b , assume i«0, x«0 
Therefore, 


— m 0 
dx dt 


di 

dt 


dL ^ dx 


SO v = Ri + L— + 


dL 2N 2 A\i 0 g 0 v 

Consider — = — —x 

dx gl - X 2 


T +dL+dx 9 2N 2 A\i 0 g 0 ■ 0 

b dX dt b (gl - X 2 )~ 


dx 

2 dt 


Thus v(t) = Ri ( t) + L-^i(t) 


taking the Laplace transform on both sides, we get 


V(s) = RI (s) + Lsl(s) 

«■ I(s) ■ Thr s ns) 


Recall equation (3) 
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Fix, i) 


s + A ^ 2l O j 


gl 


gl 


By Newton Second Law Fix, i) = m 


d 2 x 

dt 2 


„ „d 2 x A\x a N 2 ll A\lJ 1 2 I 0 
So, m = x + 7 


dt' 


gl 


gl 


Taking the Laplace Transform on both sides, we get 

g 0 go 


or 


X(s) = 

Vis) 


A\igN 2 I 0 

2 

g 0 


, A\lJJ 2 I, 
iR + Ls) ims 2 — 


• 2 T 2 
o 


gl 


( 6 ) 


A»qH 2 Iq 

glmL 


la . £ Ua , - 

L m g 0 


Let a = 


^ 2J ° b = 


glmL 


R 

T’ 


, 2 _ 


= 


m gl 


Therefore the plant dynamics of the ALPS actuator are 


» X(s) _ 

A -— 7 — r- - 


Gp{s) " V(S) is + b) is 2 ~ C 2 ) 


Which is similar to the plant dynamic equation obtained by 
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Kilgore 9 and Jayawant 3 . 

Referring to Groom 2 , the values of those parameters are, 


Io 

0.57 Amps 

A 

1. 1400918*10' 3 m 2 
— 


4 71 *10 ’ 7 H m' 1 

N 

1386 turns per coil 

m 

7.19712 kg 

g» 

0.00762 m 

R 

a 

o 

• 

00 

L 

0.1805899 h at g 0 1 


a * 20.79 
b = 44.3 
C 2 * 280.8 

Thus, the open loop transfer function is 

_ 20.79 

P s 3 + 44.3s 2 - 280.8s - 12439.59 
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d) Discussion of the plant dynamic of ASPS actuator 
The open loop transfer function is a third order, type zero, 
all poles plant system. The characteristic equation also contains 
one positive real root, so this plant is not BIBO stable. The pole 
zero diagram of the plant is shown in fig 5. In order to move the 
open loop unstable root into the stable region, we need to add 
zeros in the left half plane so that the locus are pulled into the 
stable region. In other word, a reshaping of the root locus 
(compensator) is necessary. 


pole zero locotTon of open loop plont 



figure 5 
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IQ. Compensators Design 

A. Cascade Compens ator 


The approach we used in the cascade compensator is to achieve 
the goal ( shift the root locus to the L.H.P. ) with minimum 
complexity. For the simplest case, by the knowledge that a pole 
w j_H pull the root locus to the right, a zero will pull the root 
locus to the left and a pole-zero pair close to the origin will 
decrease the steady state step error (G p is type zero) , we tried the 
general lag— lead cascade compensator 10 . 


(s 


+ 4 -) (s + 4 -) 


Gjs) = A 




( 3 + 


a r. 


) (S+-5-) 


with a =10, and gain A. 10 The lag component was fixed at 

(s+ 0 . 05 )/ (s+ 0.005) and the lead component was moved along the 
real axis. Some root locus results are shown in appendix I. After 
studying the results, we decided that we needed to increase the 
compensator complexity in order to meet the design specifications. 

Since the lag component only affected the steady state error, 
for simplicity, we tried the dual phase advanced compensator with 

a=10. 


Gjs) 


s + 
s + 


1 

T_ 

a 

T 


2 


Some root locus results were shown in the appendix II. 
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Imog Axis 




- < J ^rr> 2 


A block diagram is shown in figure 7. We selected a damping ratio 
5 = 0.7, and the maximum natural frequency. The figures of merit 


are. 


QflGHKAL PAGC IB 
OF POOH QUALITY 






poles: -86.15 ± 87 . 89 j , -93.24, -620.46 

additional gain K = 1.011406 E 6 

steady state step error = 0.06 

rise time = 0.009 sec. 

peak overshoot = 1.42 

peak time = 0.025 sec. 

settling time = 0.05 sec. ( for 5% ) 

gain margin = 16.75 dB 

phase margin = 180° 

stiffness, u n = 15 kN kg'm* 1 

stability region : 6.0939 E 4 < K < 6.8041 E 6 

The step response was shown in figure 8, and the Bode plot was 
shown in figure 9. Those results are obtained by matlab. 

figure 7. 
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step response of coscode connpensotor 



figure 8 
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Phase deg Gain dB 



figure 9 
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The 


B. Feedback Compe nsator Design 
cascade compensator will be converted to a digital 
controller later, and a computer will be involved to control the 
plant. Therefore, it is natural to design the compensator using the 
state feedback technique. This technique is flexible and convenient 
to implement. A brief derivation of the design procedure is shown 
in Appendix III. In the case where some state variables are not 
accessible, an observer (estimator) may be used. Observer design 
procedures are also shown in the Appendix III. This material are 
come D'Azzo 10 . In this design, we used the full state feedback 

technique. 

A state space representation of the open loop plant is shown 
in the figure 10. 


figure 10 

State space representation of the open loop plant 



TM 

i i 1 

it L 


O 


o 

^o.S 
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O 


I 






x^ 

4- 


J^3 



o 

o 

io.^ 


id x L ! 

O’ 


o 
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I 




The controllability matrix is, 

Mc=[b|Ab|A 2 b] 

I * 

o ! 0 • 20.79 

i i 

0 ! 20.79 ' -920.997 

i i 

20.79 ! -920.997 i 46637.99 

i i 

as det ( Me) + 0, so Me has full rank. 

The observability matrix is 
Mo=[C t |A t C t |(A t ) 2 C t ] 

1 j 0 I 0 

I I 

I I 

0 i 1 I 0 

I I 

I I 

o ! o ! i 

I I 

I I 

So Mq has full rank. Hence, this system is completely 

controllable and completely observable. 

Motivated by the performance of the cascade compensator, we 
selected the poles of the control ratio to be -86.15 + 87.89j, 

—100, which give us a good step response. A block diagram is shown 
in the figure 11. The figures of merit are, 
gain = 72.8546 E 3 
k, = 1.01 
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0.021561 


k 2 = 
k 3 = 0.0001505 

poles : -86.15 ± 87 . 89 j , -100 
steady state step error = 0 
rise time = 0.03 sec. 

settling time = 0.045 sec. (for 5 %) 
gain margin = 13 dB 
phase margin = 160° 
stiffness = 15 kNkg 1 m' 1 

A step response and Bode plot were shown in figure 12 and 
figure 13 respectivity . Those results are obtained by matlab. 

figure 11. 



29 


Amplitude 



figure 12 
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Phase deg Gain dB 




Frequency (rad/aec) 


figure 13 
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IV Computer Simulation 


We simulated the state feedback design by a four order Runge 
Kntta Method in the time domain. The Runge Kntta Method is a 
special version of the general Taylor's series expansion. 

For a general different equation. 


dX 

dt 


= J?( X ) 


Let, 


1 ! = 1 ( X ) 

r 2 = i ( x + h r x ) 

^3 = 1 ( X + h j ? 2 ) 

t \ * 1 ( X + h l 3 ) 


The next x can be approximated by 


* X * | [ ?! * 2 ( 3 

with an error of fifth power term of the Taylor series 
expansion. The h is the increment step side. 

Recall, 
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mx= 


A\1 0 N 2 (g 0 l 0 + X i) ig 0 i + l 0 X 
( 9 o "* 2 > 2 

di 


v = R i + L- 


dt 


Let 


x 1 = x = output 



v = u = input 


So 


*1 = *2 


x 2 = a 


A V ‘° N - (ff 0 Io + x 3 ) (gr Q x 3 

/n 


2 \ 2 


(?„ 2 - v> 


*> * -I * i a 


*4 - 1 

Let the control law be 

u = r - ( Jc x x x + Jc 2 x 2 + k 3 x 2 ) 


where r is the reference input. 
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need to estimate the second 


In the computer code, we 
derivative of x, which is achieved by 

(v - v„ . ) / time interval 
For the linearized plant dynamics equations, slightly 
modification of the computer codes can do the job. The results of 
computer simulation are shown in figure 14 to 19. The computer 
program was written in Pascal language, and is shown in Appendix 

IV. 


step response with linear equation 

1.2 



figure 14 


3 



figure 15 


figure 16 



figure 17 


figure 18 



figure 19 
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V. Discussion of results 


The computer simulation showed that the step response 
increased linearly with time, figure 17. That is, an unstable 
response occurred. 

Consider the plant dynamics equations, recall 


A \i a N 2 ( g 0 I 0 + x i ) ( 9 0 + J o x ) 

* = ( < 7„ 2 - x 2 ) 2 


As for | x | < | g 0 | 
1 1 


(g 0 2 - x 2 ) 2 


- 9 0 


( 1 - — 2 v 


= gr 0 - 4 [ i * <-2El> * <^) 2 -...l 2 

Jo 

= g -* [ 1 + 2 ( — ) + 3 ( — ) 2 + higher terms ] 


So 


( 9o l Q + x i) (g oi + Ip x) 


= 9a 


(g 0 2 - x 2 ) 2 
= 9 0 ~* (9 0 I 0 + x i) ( g a i + I Q x) 

[ i + 2( — ) + 3( — ) 2 + higher terms ] 
9o 2 9 2 

(g a I 2 x + g 2 I 0 i + g a x i 2 + l 0 x 2 i) 

[ l + 2(-^i) + 3 (——z) 2 + • • -1 
9o 9 0 


= J e " 4 


[ (g 0 Io 2 x ♦ g 2 Io i) 


+ (2 x 3 + 3 I a X 2 i + g 0 X i 2 ) + higher terms ] 

Jo 


Obviously, the coefficient magnitude of the third terms in the 
Taylor's series expansion is larger than the coefficient magnitude 
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of the first term. Thus, at least, we need to include the third 
terms in the compensator design. However, we cannot use the 
conventional linear design theory in this situation. 

Even thought we include the third terms in the design, the 
coefficient magnitude of the x“ terms in the series expansion 
increase without bounded ( because g c < 1 ) . So, the Taylor's series 
expansion does not exist. 
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VI. Conclusion and recommendation 

The plant dynamics equations are nonlinear, and the 
conventional linearization does not work. We recommend to design 
the compensator without linearization in the time domain. 

As, in general, if we close the loop, 

= t ( x , u ) 
at 

u = g (r , x ) 


define PI to be 

(” t (x x - r ) 2 dt 

J o 


and, minimizing the PI which is subjected to the constraint 
equation , 

^ x , g ( r , x ) ) 

dt 

, by the Langrange multiple method. The mathematics is too 
difficult to carry out, and the analysis is left for interested 
reader only. 
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InwgAxs 


G e = [ ( s + 44.3) / (s + 443) ] 2 
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Appendix ^11 

state feedback des icyn procedure 10 


Let the open-loop plant transfer function be 


g p = 


+ a 2 s‘ 


u 

+ a[ s * a 0 


X 

u 


Let, for standard notation, 


Xj. = y. 



*3 


d^y 

dt 2 


dXj. _ d£ 
" dt dt 


2 ' 


dx 2 _ d}y _ dx^ = d 3 y 

dt ~ dt 2 3 ' dt dt 3 


and 


y = 


dX 

dt 

[ 1 


= A X + h u 
0 0 ] ^ = £ r 2C 


Let the control law be 


u = r - k T X 

where k T * [ -fc 2 ^3 1 


Consider 
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H(s) 


k T X(s) _ k T X 
X(s) c 7 X 


H(s) 


x x + k 2 x 2 + Jc 3 x 3 

x i 

K X 1 + k 2 SX 1 + *3 S 2 *! 
*1 

= Jc x + sie 2 + s 2 Jc 3 


C g (*! + gJc 2 S 2 Jc 3 ) 
s 3 + a 2 s 2 +a 3 s +a 0 


and the overall transfer function is 

Y(s) _ G _ C ° 

R{s) ~ 1 + G H (s 3 + a 2 s 2 + a x s + a 0 ) + c 0 (ic 3 s 3 + k 2 s + k x ) 

c o 

s 3 + (a 2 + c 0 Jc 3 ) s 2 + (a 3 + c Q k 2 ) s + (a 0 + Co^) 


using the final value theorem, 


lim y ( t) : = y aa = lim s Y(s) 

t — 0 -o 

3 Cg RjS) 

y sB < c) = S 3 7 (a 2 + c 0 Jc 3 )s 2 + (a 2 + c 0 k 2 ) + (a 0 + c 0 ^) 


for step input R(s) ■ s' 1 
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••• y SB = 


a 0 + C 0 J c,_ 


: = 1 for zero error 


ie, k, = 1 - 


which is fixed. By appropriate selecting the value of k T , we can 
implement any desired characteristic equation as we want. 


original plant: 


observer design 


dX 

dt 


= A X + h u 


y = Z T X 


Let k be the estimated state vector. 
Let 
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-p = AX + hu+L { y - 9 ) 

at 

9 = c 7 X 

where L is the observer matrix 

L = l 1 i i 2 J 3 ) 

define X - X 

ds. dX _ d ^ 

~*~dt ~ dt dt 

= A X + h u - (A2 + hu+L( c T X - c. T 2) ) 

= A( x - 2) - L £ r ( x - 2 ) 

= ( A ~ L Q t ) £ 

By appropriate selecting the eigenvalues of e f the error of the 
estimated state vector will died out very quick. A state diagram is 
shown in figure 20. 


figure 20. 













Appendix IV 



{$N+,E+} 

program project (outputo) ; {3\21\1993} 

var 

d : char ; 

i, j, n, nstep, e, b: integer; 

h, h2 , Xlmax, X4max, t, r, u, dum:real; 

k:array[l. .3] of real; 

x,y:array[l. .2,1. .4] of real; 

f : array [1. .4,1. .4] of real; 

outputo : text ; 

{ This is a Runge-Kntta method of order 4 . > 
procedure initizing; 


begin 


t:=0; 

{initial time} 

h:=l/1024; 

{incremental time step} 

n:=4 ; 

{number of equations} 

for i:=l to 2 do 

{initial value} 


for jj^l to 4 do 
begin 
x[i, j] :=0; 

y[i, j] := 0 ; 

end; 

for i:=l to 4 do 
for j:=l to 4 do 
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nstep: =4096 ; {number of step} 

r :=1; 

k[l] :=1. 01; 
k[2] :=0. 021561; 
k[3] :=0. 0001505; 
end; 

procedure get_value_f; 
begin 

f [e, 1] := x [ b , 2 ] ; 

f [e , 2 ] := (382 . 3996e-6) * (4 . 3434e-3 + 

x[b,l]*x[b,3])*(0. 00762*x[b, 3 ] + 

0 . 57 *x[b, 1 ] )/ (sqr(sqr(0. 00762) - sqr (x[b, 1] ) ) ) ; 
dum := (y[b, 2] - x[b,2])/h; {estimated the second derivate} 
u ; = r - (k[ 1] *x[b, 1] + k[2]*x[b,2] + k[3]*dum); 

{control law} 

f [e,3] := -44 . 2995*x[b, 3 ] + 5.5374*u; 

f [e,4] := 1; { x[4]= time } 

end; 

procedure RK4SYS; 
begin 

h2 :=0. 5*h; 

for j:=0 to nstep do 
begin 
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for e:=l to 2 do 


begin 

b:=e; 

get_value_f ; 
for i:=l to n do 

x[2,i]:=x[l,i] + h2*f[e,i]; 

end; 
e:=3 ; 

get_value_f; {get f3} 
for i:=l to n do 

x[2, i] :=x[l, i] +h*f [e, i] ; 

e : =4 ; 

get_va lue_f ; { get f 4 > 

y[l,2]:-x[l,2] ; 

Y [ 2 , 2 ] : =x[ 2 , 2 ] ; 

for i:=l to n-1 do {compute next x(t+h)} 

x[l, i] :=x[l / i]+h*(f[l / i]+2*(f[2,i]+f[3,i] )+f [4, i] ) /6; 

x[l,4]:= t+j*h; {advance solution} 

if j mod 64 =0 then {write the result} 

begin 

for i:=l to n do 
write ( outputo , x [ 1 , i ] , ' , ') ; 
writeln(outputo) ; 
writeln(x[l, 1] ,','/X[l,4]); 
end {if loop} 
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end; {for j loop} 
end; {RK4SYS} 

begin {main} 

assign (outputo, 'a:\pl.dat' ) ; 

rewrite (outputo) ; 

writeln (outputo) ; 

initizing; 

repeat 

RK4SYS; 

write('want change Y/N ?'); 
readln(d) ; 
if d='y' then 
begin 

writeln ( 'kl=' ,k[l] , ' k2=',k[2],' k3=',k[3]); 
write ( ' enter kl,k2,k3 / ); 
readln(k[l] ,k[2] ,k[3]) ; 
end; 

until (do'y'); 

writeln (outputo, ' Job completed. '); 
writeln ( ' Job completed . ' ) ; 
close (outputo) ; 
end. 
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